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In this review we summarise recent results for the complex eigenvalues 
and singular values of finite products of finite size random matrices, their 
correlation functions and asymptotic limits. The matrices in the product 
are taken from ensembles of independent real, complex, or quaternionic 
Ginibre matrices, or truncated unitary matrices. Additional mixing within 
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and some straightforward generalisations. The asymptotic results for large 
matrix size include new microscopic universality classes at the origin and 
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other parts of the spectrum the known standard universality classes have 
been identified. In the limit of infinite products the Lyapunov and stability 
exponents share the same normal distribution. To leading order they both 
follow a permanental point processes. Our focus is on presenting recent 
developments in this rapidly evolving area of research. 
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1. Introduction 

The study of products of random matrices goes back to the early days of ran¬ 
dom matrix theory, when in 1960 Furstenberg and Kesten [1] studied them 
in the context of dynamical systems and their Lyapunov exponents. When 
multiplying matrices generically the product becomes non-Hermitian. Hence 
it is natural to choose each factor independently to be non-Hermitian, taken 
from a Gaussian distribution in the simplest case, the Ginibre ensembles [2] 
with real, complex or quaternionic matrix elements. These are the analogues 
of the three classical Wigner-Dyson ensembles, labelled by /3 = 1, 2,4 respec¬ 
tively. When modelling unitary time evolution the choice of matrices to be 
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multiplied is slightly less trivial. Taking factors from ensembles of Haar dis¬ 
tributed orthogonal, unitary or symplectic matrices does not lead out of the 
one matrix case, due to the invariance of the Haar measure. For the product 
to become non-trivial instead one may choose to multiply truncated unitary 
matrices which are distributed according to the Jacobi measure, leading to 
a sub-unitary evolution. Other choices are of course possible, e.g. by multi¬ 
plying unitary matrices that are not Haar distributed [3, 4, 5]. 

In any case one has two choices in studying spectral properties of the 
product matrix: either the eigenvalues which are generally complex, or the 
singular values which are real positive. The latter were studied hrst in order 
to dehne and study the Lyapunov exponents. Their application to dynamical 
systems has lead to the activities summarised in [6]. Examples for more 
recent applications include wireless communication [7] and combinatorics [8]. 
Although in general complex eigenvalues and singular values are not related 
individually, in the limit of inhnite products the radii and singular values 
become identical, in other words the stability and Lyapunov exponents have 
the same normal distribution [9, 10], as hrst conjectured in [11]. 

It is quite surprising that the spectral properties of products of random 
matrices were hrst determined in the inhnite product case, using mainly 
probabilistic tools, and in the limit of inhnite matrix dimension (keeping the 
number of factors hxed) using planar Feynman diagrams [12] and probability 
theory [7, 13, 14], see [15] for a very recent review on the free probability 
approach. The common point to the latter two approaches is that they 
contain information about global spectral properties. These were found to 
be universal when multiplying matrices from various different ensembles [16], 
including analytical results for factors from the elliptic ensemble [17]. 

Apart from special cases for 2x2 matrices, see e.g. [18], explicit results 
for the product of M matrices of size N x N with M and N hnite are 
very recent, starting from [19, 20, 21] for the complex eigenvalues and from 
[22, 23] for the singular values. They reveal a determinantal and Pfafhan 
structure and will be the main topic of this review. Inherited from the 
fact that the product of normal or gamma random variables is distributed 
according to the Meijer G-function [24], these special functions appear in the 
weight and determinantal expressions for the products of random matrices 
too. Obviously the detailed knowledge of all eigenvalue and singular value 
correlation functions has opened up the possibility to study local spectral 
properties, and we will also review the recent progress on hnding known and 
new universality classes. 

In [19] the joint density of complex eigenvalues for products of M ma¬ 
trices of size N X N was derived and the kernel of orthogonal polynomials 
(OP) in the complex plane that determines ah correlation functions was 
given. While in the bulk and edge large-A^ scaling limit at hxed M the 
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respective Ginibre universality classes for a single matrix were recovered, at 
the origin a new universality class labelled by M with a hypergeometric ker¬ 
nel appeared [19]. These results were generalised to products of rectangular 
matrices (these are also called induced Ginibre matrices [25]) for /3 = 4 in 
[21], to the rectangular /3 = 2 case [26], and to the joint density for (3 = 1 
in [27]. Here also products of rectangular truncated unitary matrices and 
products of mixed type were considered, for which a weak commutation re¬ 
lation was proved. In all these results the change of variables from matrices 
to complex eigenvalues uses a generalised Schur decomposition. 

This decomposition was well-known for M = 2 and /3 = 1, 2, see e.g. [28], 
whereas its generalisation to an arbitrary number of matrices M was given 
in [29, 26] for (3 = 2, and was extended to /3 = 4 in [30]. A similar ex¬ 
tension exists for (3 = 1 and was implicitly used in [31, 27]. In [26] also 
mixed products of Ginibre and inverse Ginibre matrices, and of truncated 
unitary matrices and their inverses were considered, and their joint densities 
and corresponding kernels were given. The large-limit for such products 
of unitary matrices truncated form U{N + k) to U{N), also called random 
contractions, was studied in [32]. At strong non-unitarity in the bulk and 
edge scaling limit the universal result for a single Ginibre matrix was recov¬ 
ered. In the origin limit the same new class as for products of M Ginibre 
matrices [19] mentioned above was found. Very recently, a more rigorous 
derivation of this bulk and edge universality was presented in [33], including 
the case of products of rectangular Ginibre matrices when both the inner 
and outer edge become soft. At weak non-unitarity the result for a single 
matrix [34] was extended to a new kernel labelled by Mk in [32]. A discus¬ 
sion for M > 1 products of truncated orthogonal (and unitary symplectic) 
matrices was given in [27], albeit an understanding of the universal kernels 
is still lacking, for a single truncated orthogonal matrix see however [35]. 

The distribution of radii of complex eigenvalues of products of Ginibre 
matrices is given by permanents [20, 30], cf. [36] for M = 1. The corre¬ 
sponding gap and overcrowding probabilities were derived in [20] for (3 = 2 
and in [30] for /3 = 4, including the corresponding asymptotic expansions 
for finite and infinite point processes, for fixed M. The limiting distribution 
of the largest radius of all eigenvalues was shown to interpolate between the 
log-normal and Gumbel distribution (valid for a single Ginibre matrix [37]) 
in a double scaling limit in [38], where both N and M become infinite. 

Turning to singular values their joint density and corresponding kernel of 
orthogonal functions were derived in [22] for square and [23] for rectangular 
matrices with (3 = 2. These enjoy a relation to multiple OP [39] and a 
universal correlation kernel labelled by M was found in the local origin 
scaling limit [40]. For M = 2 it agrees with the limiting kernel in a Gauchy 
two-matrix model [41], a correspondence that was extended to the Gauchy 
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multi-matrix model and general M very recently [42], cf. [43]. This result 
was generalised in [44] to a mixed product of Ginibre matrices times a single 
truncated unitary matrix, introducing the notion of polynomial ensembles. 
There, their relation to a certain type of biorthogonal ensembles previously 
studied in [45, 46] was pointed out. Polynomial ensembles enjoy special 
invariance properties, related to products of random matrices [47]. The 
same limiting kernel as in [40] was also found in [48] for products of Ginibre 
and inverse Ginibre matrices, where the average characteristic polynomial 
was computed as well. The zeros and asymptotic analysis of the average 
characteristic polynomial was performed in [49] . The Fredholm determinant 
for the gap probability at the origin was shown to satisfy a system of non¬ 
linear ordinary differential equations [50] . In the bulk and at the soft edge it 
was conjectured in [22] to find the universal sine- and Airy kernel respectively, 
see [51] for a very recent proof including mixed products of Ginibre matrices 
and their inverse. The correlation functions of singular values of truncated 
unitary matrices is also currently under way [52]. So far information on local 
properties of the singular values for the /? = 1,4 classes has been very difficult 
to access, due to absence of the corresponding Harish-Ghandra integral. Very 
recent progress has been possible using supersymmetric techniques, see [53]. 

The limiting positions for the Lyapunov exponents of products of /3 = 2 
[54] and /3 = 4 [55] Ginibre matrices were derived using probabilistic meth¬ 
ods, including correlated Gaussian distributions for each factor. The detailed 
knowledge for the joint densities described above allowed to show [32] for 
uncorrelated complex Ginibre matrices, that each Lyapunov exponent be¬ 
comes normally distributed, following a permanental point process. The 
same leading order behaviour holds taking the identical limit of the radii of 
complex eigenvalues, called stability exponents. This leads to a one-to-one 
correspondence between limiting radii and singular values. The expressions 
for the stability exponents were extended to /3 = 4 and /3 = 1 most recently 
in [10], where the latter case was obtained under the assumption that for 
large M all eigenvalues of the real product matrix become real. Such a 
behaviour was previously observed in [56] and proved for general N in [31]. 

All aforementioned finite-V and -M results have assumed that each fac¬ 
tor in the product is independent^. To date only complex eigenvalues of 
products of M = 2 matrices that are coupled through an Itzykson-Zuber 
term and thus not independent have been studied for finite-V at /? = 2 [57], 
/3 = 4 [58], and /3 = 1 [59]. They are given in terms of parameter depen¬ 
dent families of Laguerre polynomials in the complex plane, and we refer to 
the review [60] where all these results are summarised. Results for singular 
values of products of such dependent matrices are currently under way [61]. 


^ Of course powers of a single matrix form an exception, leading to the same global 
density as the product of independent matrices, see e.g. [13]. 
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The remaining content of this paper is organised as follows. In section 2 
we will present exact results for finite N and M including the joint densities, 
kernels and correlation functions; the section is divided in to two subsections 
2.1 and 2.2 discussing results for complex eigenvalues and singular values, 
respectively. In both cases, we consider Ginibre and inverse Ginibre matrices 
as well as truncated unitary matrices. We turn to the local large-limits 
at fixed M depending on the location in the spectrum in section 3. Here 
we will be very brief and only give more details for the limiting kernels that 
were not previously known. Again, we will hrst discuss complex eigenvalues 
in subsection 3.1 and then discuss singular values in subsection 3.2. The 
following section 4 is devoted to the opposite infinite product limit at fixed 
N. Our discussion of open problems follows in section 5. 

2. Exact results for finite products of finite size matrices 

2.1. Complex eigenvalues 

2.1.1. Ginibre matrices 

We begin with the simplest case: A product of rectangular Ginibre matrices. 
We emphasise that this special case is very illustrative, since the main ideas 
from the treatment of products of Ginibre matrices extend to all the other 
examples described in this review. 

We are seeking the statistical properties of complex eigenvalues of the 
product matrix 

Hm = XmXm-i-■ ■ X 2 X 1 ( 2 . 1 ) 

where each Xj is an Nj x Nj-i matrix distributed according to a Gaussian 
density 


n \ gNjNj-i/2 


Tr wwj 
27 ^ f 


j = l,...,M. 


( 2 . 2 ) 


Here the matrices Xj [j = 1,..., M) are independent and the index fd = 
1, 2,4 denotes whether real, complex or quaternionic matrices are considered 
(we only multiply matrices with the same /3). The parameter 


7 = 



for 13 = 1,2 
for f3 = A 


(2.3) 


is related to the fact that quaternions are given by their 2x2 matrix rep¬ 
resentation, which implies that the eigenvalues come in complex conjugate 
pairs (see e.g. [62]). If we disregard eigenvalues which are trivially zero, then 
we can choose N = Nq < Ni < ■ ■ ■ < Nm without loss of generality; this is 
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due to a weak commutation relation [27]. Furthermore, we can parametrise 
the matrix product as 


XmXm-1 • • • X 2 X 1 = U (^^mXm-1^- • • ^ 2 X 1 ^ ^ ^2.4) 

where U is an Nm X Nm orthogonal (/3 = 1), unitary (/3 = 2), or unitary 
symplectic (/3 = 4) matrix and each Xj is an N x N matrix. This parame- 
terisation results in the measure [27] 

M M 

n dX,P{Xj) = dfi{U) n dXjP,^ (Xj) , (2.5) 

i=i i=i 


where dXj and dXj are the flat measures over all independent matrix ele¬ 
ments, d^{U) is the normalised Haar measure, and each Xj is distributed 
according to the induced density 


P,,{Xj) (X det[X,xj]^"^/(2T') exp 


27 


TtXjX] 


( 2 . 6 ) 


Here we have introduced the convenient notation Vj = Nj — N for the dif¬ 
ferences between matrix dimensions. We emphasise that (2.5) enables us 
to reduce the problem involving rectangular matrices to the problem involv¬ 
ing square matrices, which is a considerable simplification (we refer to [27] 
for a more thorough discussion). It should be noted that this simplifica¬ 
tion is possible for all examples mentioned in this review. Finally, we note 
that the induced densities (2.6) are isotropic, i.e. they are invariant under 
Xj — >■ UXjV, where U and V are orthogonal (j3 = 1), unitary (/3 = 2), 
or unitary symplectic (/3 = 4) matrices. This is an extremely important 
observation needed for the weak communication relation [27], which states 
that any averaged property of such a product matrix (depending only on 
the product matrix itself, not on the individual matrices) is independent of 
the ordering of the factors. For products of rectangular Ginibre matrices 
this implies that all averaged properties (e.g. correlations between eigenval¬ 
ues or singular values) are invariant under permutations of the indices 1/^ 
(m = 1,..., M). The reader might note that all products we consider in this 
review are independent of the ordering of the factors. This originates from 
the fact that all ensembles we consider are isotropic. 


* * * 


We are now ready to state the normalised joint probability density function 
(jpdf) of complex eigenvalues for the product matrix H^- For (3 = 2 the 




jpdf reads [19, 26, 27] 


n (2-7) 

^Ar,i/ n=l l<j<Z<Ar 


Apart from the corresponding weight function given by a Meijer G-function 


= G 


M,0 

0,M 




’ 2 





( 2 . 8 ) 


it agrees with the jpdf for a single Ginibre matrix M = 1, where the repulsion 
of eigenvalues is given by the squared absolute value of the Vandermonde 
determinant. The normalisation constant is given for all three /3 = 1, 2,4 by 


^N,u = 2(2-/3)7MAf(Af+l)/4 11 11 ^ ( 2 

n=l j=l ^ 

The following steps are needed to derive the jpdf in (2.7) for /? = 2 for general 
M (and subsequently for the other /3’s): i) a generalised Schur decomposition 
[19, 29, 26] of the individual factors distributed according to (2.6), and ii) 
the fact that the upper triangular matrices from this decomposition decouple. 
The remaining integrals can be performed [19, 26, 27] and lead to the Meijer 
G-function as the weight in (2.8). For the most general set of indices the 
Meijer G-function is defined as [63] 

\ ^ f ds njLi - s) ll]=i r(l - Oj + a) 

7 “ 27rz Jc irj=m+i r(l - + s) h^n+i r(ai “ «)' 

( 2 . 10 ) 

The integration contour C depends on the poles of the Gamma functions, cf. 
[63]. In the simplest case with p + q < 2{m + n) and z in the upper half 
plane it runs from —zoo to +zoo leaving all poles of T{bj — s) to the right 
and of r(l — Oj• + s) to the left. Empty products are defined as unity. For 
the special case (2.8) an alternative multiple integral representation of the 
Meijer G-function exists [19, 26] 


^m,n 


Cii , • • • 5 dp 
bi,...,bq 



( 2 . 11 ) 


Note that the right hand side corresponds to the density of a product of 
(unnormalised) scalar-valued random variables. Here we get the explanation 
why the Meijer G-function appears: it was previously known for products 
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of Gaussian or gamma random variables that their products is distributed 
according to the Meijer G-function [24], This property is passed on to prod¬ 
ucts of random matrices as well. It follows immediately from the represen¬ 
tation (2.11) that the (bi-)moments of the weight functions are given by 


fz w^~'^ 


{z) z 


k 


M 

= JJr( 

i=i 


r'j k 1). 


( 2 . 12 ) 


Here the moments are zero for k ^ i since the weight function (2.8) is 
invariant under rotation in the complex plane. 


* * * 


The k-point density correlation functions are defined following [62], 

Rk{zi,...,Zk) = d^Zk+i ■ ■ ■ d^ZNV^p^fizi, zn). (2.13) 

We will consider jpdfs with the following form: 

1 ^ 

= (2.14) 

n=l l<£<k<N 

where wl^^‘^{\z\) is an arbitrary weight function that only depends on the 
modulus. Then it is well-known that the corresponding OPs are monomials, 


d^zw'^ '^{\z\)z’'z*^ = 6k,ihk. 


(2.15) 


Furthermore, it follows from general considerations that the correlation func¬ 
tions including the jpdf for k = N form a determinantal point process [62] 


^{zi,...,zk) = Ylw^ ^ {\zj\) K^^{zi,Zj) 

j=i 


where the kernel is given by 




N-1 


{zu*y 


(2.16) 


(2.17) 


n=0 


Here denotes the squared norm determined through the orthogonality 
relation (2.15). The normalisation constant in (2.14) is uniquely determined 
by these squared norms as well. Explicitly, we have 


N-l 


n 


n=0 


(2.18) 
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For the product of complex (^0 = 2) rectangular Ginibre matrices the jpdf 
is given by (2.7) with weight (2.8) and the squared norms are immediately 
obtained from (2.12). Thus the kernel for the product of rectangular Ginibre 
matrices is given by [19, 26] 


N-l 




{zu*Y 


n=0 0 ™=! + 1 + l^m) 

while the corresponding normalisation constant is 

M 


(2.19) 


4 ,. = Nl7r 


M 


nr(^ + i + ^^') • 

i=i 


( 2 . 20 ) 


This determines all A:-point correlation functions for (3 = 2 via (2.16) with 
weight (2.8). 

The second set of results that follows from eq. (2.14) for /3 = 2 is the 
joint density of the radii rj of the complex eigenvalues, zj = rje^'^G It is 
obtained by integrating over all angles 6j leading to a permanent [20]: 


N 


n 

n=l 





n"i2>ru.'>-=(r„) 


2 


per 


„2l-l 


( 2 . 21 ) 


where we have included the factors from the radial measure dr^rn- The 
radii thus become independent random variables, generalising the result in 
[64]. Moreover, the hole probability that a disc of radius r centred at the 
origin is empty of eigenvalues becomes 


Prob[Vj : rj > r] 



dOnV^^^^{zi,...,ZN) 


N 


n 

n=l 


^M+1,0 


( 1 

nf=inn+i^j) 



( 2 . 22 ) 


(2.23) 


The second line was derived in [30]. It agrees with [62] for M = 1 with 
1^1 = 0. An alternative result for M = 2 was previously obtained in [65]. 


* * * 


The quaternionic k-point correlation functions for ,0 = 4 are dehned identi- 
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cally to the complex case (2.13). Here we will consider jpdfs of the form: 

1 ^ 

= —^Ww^-^{\Zn\)\Zn-zl\^ JJ \zi-Zk\^\zi-zl\^, 

n=l l<k<^<N 

(2.24) 

where the weight function is invariant under rotation in the complex plane. 
We restrict ourselves to the upper half plane C+, due to complex conjugated 
pairing of eigenvalues for /3 = 4. It should be noted that except for the 
weight function the jpdf is identical to that of a single Ginibre matrix M = 1. 
Additionally, we define the skew-symmetric product 

{f,9)s^ [ d^zw^^‘^{\z\){z* - z){f{z)g{z*) - f{z*)g{z)). (2.25) 

Jc+ 

It follows from general considerations that the corresponding A:-point corre¬ 
lation functions including the jpdf form a Pfaffian point process [62, 66] 




n 

1=1 


w^=\zj)iz*-z,) Pf 


Kl^^-\zi,z*) -K^^-\z*,z*) 


K^'^{zi,Zj) -K^^{z*,Zj) 


(2.26) 


with the kernel given in terms of monic skew-orthogonal polynomials 


r^f 3 = 4 r ^ ^ P2n+l{z)p2n{u) - P2n+l{v)p2n{z) 

^N.u = 2^ - - -. (2.27) 

n=0 


Here the skew-orthogonal polynomials Pn{z) are defined such that they satisfy 
the skew-orthogonality relations [62, 66] 

{P2k+l,P2i+l)s = {P2k,P2i)s = 0 ) {P2k+l,P2i) s = (2.28) 

We know that the weight function in the jpdf (2.24) is invariant under rota¬ 
tion in the complex plane. This implies that 



dh w^=\\z\)z’^z 


*e 


^k ^k,i 1 


(2.29) 


where Sk are constants depending on the weight. It was pointed out in [21] 
that in order to find the skew-orthogonal polynomials it is sufficient to cal¬ 
culate the constants, Sk- Explicitly, we have 


P 2 n{z) 


n 


E 


H 


'-i=k+i 


S2l 


^2k 


2n-\-l 


S2e-1 


P 2 n-i-l{z) = Z‘ 


hn = 2s2n+l- (2.30) 
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The main idea behind this result is that the odd polynomials are always 
monomials, due to the rotational invariance of the weight. 

We can now turn to products of rectangular Ginibre matrices with quater- 
nionic matrix elements (/3 = 4), where we refer to [21, 27, 30] for details. 
The derivation of the jpdf for the complex eigenvalues is based on the gen¬ 
eralised Schur decomposition for quaternionic matrices, see [30]. With the 
notation given above, we have 

1 ^ 

JJ \zi - Zk\‘^\zk-Z* , 
^N,u n=l l<fc<Z<Ar 

(2.31) 

with weight and normalisation constant defined in (2.8) and (2.9), respec¬ 
tively. Inserting (2.8) with /3 = 4 into (2.29) yields 

M 

= 2 Min+i)+i n + n + l). (2.32) 

m =1 


This determines the skew-orthogonal polynomials (2.30) and therefore all cor¬ 
relation functions via (2.26), with weight and kernel given by (2.8) and (2.27), 
respectively. 

Turning to the distribution of radii we obtain up to factors of 2 the same 
result as in eq. (2.21) after integration 



• • • ) ^n) 


n 


N 

n=l 


AttwH, ‘^{rn) 


7/3=4 


per 


„4i-3 

^i-1 


(2.33) 

This is true despite the initial repulsion of eigenvalues from the real line in eq. 
(2.31). The radii are once again independent random variables, generalising 
the results of [37]. Likewise the /3 = 4 hole probability defined as in eq. (2.22) 
was obtained in [30], 


1 

N o,2n-F2i/i,..., 

Prob[Vj : J'j > "r] = ]^-^ 


‘Ifl -|- ‘ll^M 


2M^2 


Y\f=in‘^n + 2v, 


n=l 


(2.34) 

For M = 1 it agrees with [62]; for M = 2 see [65] for a different expression. 


* * * 


Finally, we state the jpdf for products of rectangular real Ginibre matrices 
{13= 1). It is well-known that the eigenvalues of a real matrix are either real 
or come in complex conjugate pairs. The main difficulty for real matrices is 
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that a complete triangularisation is possible if and only if all the eigenvalues 
are real [28] . Typically, a real Ginibre matrix will have both real and complex 
eigenvalues which prevents this. However, it is always possible to make an 
incomplete triangularisation involving 2x2 matrices. If N is even then it is 
possible to write down the jpdf in terms of 2 x 2 matrices [27], 


N/2 

iPfp=J,(Zi,..., Z^/2) = ^ n n |det[Z,®l2-Il2®^fc]|, 

N,v n=l i<k<l<^ 


(2.35) 

where each Zn is a real 2x2 matrix, hence its eigenvalues are either real 
or a complex conjugate pair. These eigenvalues can be identified with the 
eigenvalues of the original N x N matrix H^- A similar expression for N 
odd was given in [27], but will not be repeated here. The weight function 
with a matrix argument is given by 



det Zj I e 


-i Zi), 


(2.36) 


which should be compared to the expression (2.11). The fact that the 
jpdf (2.35) is only known up to 2 x 2 matrix integrals plus the fact that 
one has to distinguish real and complex conjugate eigenvalue pairs has so 
far prevented progress in computing density correlation functions or the dis¬ 
tribution of radii. In [27], a general approach was presented using that the 
eigenvalues of a 2 x 2 matrix can be linked to its singular values. This 
approach led to an M-fold integral representation of the two-point weight, 
but the expression was too complicated for any further calculations to be 
tractable. 

As mentioned above, a complete triangularisation is possible when all 
eigenvalues are real, Zj = Xj G M (j = 1,... ,N). In this special case, the 
approach used for /3 = 2,4 can be extended in a straightforward manner 
leading to a jpdf with all eigenvalues real [31], 




(xi,.. .,xn) 


Z'N,y n=l l<j<l<N 


(2.37) 


with the weight and normalisation given in eqs. (2.8) and (2.9), respectively. 
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If N is even then the probability that all eigenvalues are real is given by [31] 


N 

n 

n=l ■ 


1 ’ 

det 

i<fc<£<l 


'-^M+1,M+1 I 


/3 _ « 

12 2 

EL 
’ 2 


3 _ I'M 
2 ''•:••• 1 2 2 
0, ^ + /u, + k 


- ^ - ^,1 


nlin“=ir(^ + f) 

A similar expression holds for N odd. 


(2.38) 


2.1.2. Ginibre and inverse Ginibre matrices 

Next we consider the generalised eigenvalue problem given by the character¬ 
istic equation 

det[Yi ■■■YlX-Xm ■■■Xi] = 0 , (2.39) 

where each Xj is an Nj X Nj-i matrix and each Yj is an Nm+j-i x Xm+j 
matrix. If each Xj and Yj is distributed independently according to (2.2), 
then the generalised eigenvalue problem may formally be thought of finding 
the eigenvalues of a mixed product of Ginibre and inverse Ginibre matrices, 

{Yi---Yl)-^Xm---Xi. (2.40) 


This was studied for square matrices in [26] ; the jpdf of complex eigenvalues 
reads 


.p/3-2 

jpdf,I/,/i 


{zi, ...,zn) 


n=l 


with weight function 


w 


0=2 

U,fl 


(^) 


^M,L f-N - ^ 1 , ... , -N - fiL 



(2.41) 


(2.42) 


Here we use the notation Vi = Ni — N for i = 1,... ,M and pLj = NM+j — X 
for j = 1,... ,L. Equivalently the weight function can be written in the 
following integral representation 




= vr 


n 

l=l ■ 


d^ U£e 




\U£\ 




M 


xH 

m=l 


d^Zr, 


p— 


2^*771 ^2 



ZM---ZI 

UL-'-Ui 


(2.43) 
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It follows that the squared norms defined via (2.15) are given by 

M L 

= vr n + n + 1) n r(iV + - n). (2.44) 

m=l i=l 


Thus the normalisation constant directly follows from (2.18), while all k- 
point correlation functions are given via (2.16) with weight (2.42) and ker¬ 
nel (2.17). Note that in both (2.18) and (2.17) we have 0 < n < 1, which 

ensures that the gamma functions in (2.44) are well-defined. For L = 0, the 
jpdf (2.41) reduces to the previous results (2.7), while for M = L = 1 it 
reduces to the spherical ensemble introduced and solved in [67]. 

The joint density of the radii is given by (2.21) and a simple calculation 
starting from the definition (2.22) yields the hole probability 


N 

Prob[Vj : ’’j > ^] = ]^ 

n=l 


*^L+1,M+1 


I'n - N - m,... ,n - N - HL,! 

Y 0,n + ni,... ,n + VM 



n^i r(n + + 1) nLi r(iV + fMe-n) 


(2.45) 


* * * 

The result for (3 = 2 can be extended to /3 = 4 in a simple manner. The jpdf 
of complex eigenvalues is of the same form as (2.24) except that the weight 
function and the normalisation constant changes. The weight function for 
/3 = 4 can be expressed in terms of (2.42), 

(2.46) 

It follows from (2.29) that 

M L 

n r(2z^m + n -|- 1) r(2A^ -|- 2fijn — n -|- 1). (2.47) 

m=l £=1 

This determines the skew-orthogonal polynomials (2.30) and therefore all 
correlation functions via (2.26) with the kernel given by (2.27). The corre¬ 
sponding spherical ensemble (M = L = 1) for /3 = 4 were studied in [68]. 
For /3 = 1, the spherical ensemble was studied in [69] but we will not discuss 
the generalisation to arbitrary M and L for /3 = 1 here. 

2.1.3. Truncated unitary matrices 

Now we turn to products of truncated unitary matrices. Consider M inde¬ 
pendent orthogonal (/3 = 1), unitary (/3 = 2), or unitary symplectic (/3 = 4) 
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matrices Uj of different sizes Kj for j = 1,... , M. Let the unitary matrices 
be uniformly distributed with respect their corresponding Haar measure; we 
seek the truncation of each matrix Uj to its upper-left sub-block Xj of size 
Nj X Nj-i. We are interested in the complex eigenvalues of the following 
product of M truncated unitary matrices 

\Im^Xm---Xi. (2.48) 

For Kj — Nj — Nj-i > 0, projection of the Haar measure on the original 
group to a measure for the sub-block gives 

Pj{Xj) oc 0[1 - x]Xj] det[l - ^ (2.49) 

for each Xj. Here 0 denotes the Heaviside theta function of matrix argument. 
In the more general case see [32] for an integral representation of the measure. 

For /3 = 2, the jpdf of complex eigenvalues derived in [26, 27, 32] is of 
the same form as eq. (2.7), with the weight function given by 


w 


/3=2 

1/,K 


^M,0 

'^M,M 


fni,... ,Km 



(2.50) 


Here we have Vj = Nj — N and Kj = Kj — Nj-i, with the restriction 
Kj — Vj > 0 for j = 1,..., M. The theta function in (2.49) is included in the 
Meijer G-function, because the weight function (2.50) is strictly zero outside 
the unit disc. For example we have [63]: 


^1,0/I 

‘-mIo 



0(1-|zp) 


and G 


2,0 

2,2 


1,1 

0,0 


= — log] 2 ;|^ 0 (l — \z\^). 


(2.51) 

Again the weight function can be written via an M-fold integral representa¬ 
tion [26, 27, 32], 


M 







d^Zr, 


I 21^77 


(1 


— 1 
I ) 


{z-ZM---Zi), (2.52) 


where the integration domain, D, is the unit disk. From this representation 
we hnd the squared norms 


h 


n 


M 


=^n 

m=l 


r(^'rn + + 1) 

T{Km + n+l)' 


(2.53) 
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The kernel and the normalisation constant are uniquely determined by (2.17) 
and (2.18), respectively. 

In complete analogy to the Ginibre case we may go to the generalised 
eigenvalue problem (2.39) and consider a mixed product of truncated and 
inverse truncated matrices. We introduce L additional independent unitary 
matrices Vj of different sizes Tj for j = 1,...,L, which are distributed 
uniformly with respect to Haar measure. Truncating each matrix Vj to 
its upper-left sub-block Yj of size Nm+j-i x Nm+j, the general eigenvalue 
problem gives rise to a jpdf with exactly the same structure as (2.7), except 
that the weight function (2.50) is replaced by [26] 

/3=2 f \ _ r<^,L I |2^ 

^ _]V — n, . . . , — j 

(2.54) 

Here pj = NM+j — dY and Tj = Tj—Nu+j-i satisfy the restriction Tj—pj > 0 
for j = 1,... ,L. Note that if M = 0 then the weight function (2.54) is 
strictly zero inside the unit disc. 


* * * 

The result for /3 = 2 can be extended to /3 = 4 in a simple manner. The 
jpdf of complex eigenvalues is of the same form as (2.31) except that the 
weight function and the normalisation constant changes, see [27]. The weight 
function for the product of truncated unitary symplectic matrices can be 
expressed in terms of (2.50) as 


w 




It follows from (2.29) that 




^ -r-r r(2r'm -|- u -|- 1) 

^-l-i r(2Km + n) 

m=l ^ ' 


(2.55) 


(2.56) 


This determines the skew-orthogonal polynomials (2.30) and therefore all 
correlation functions via (2.26), with the kernel given by (2.27). 


2.1.4. Mixed products 

Finally, we mention that no new techniques are required in order to study 
more complicated products constructed from different combinations of Gini¬ 
bre, inverse Ginibre, truncated unitary as well as inverse truncated unitary 
matrices. Such products will have a jpdf with a structure similar to the 
examples above and the weight function can again be expressed as a Meijer 
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G-function, albeit with more indices. In fact, products of Ginibre and trun¬ 
cated unitary matrices have previously been studied in [27] for all /3, although 
the most general results are restricted to /3 = 2,4 due to the incompleteness 
of a generalised real Schur decomposition as discussed in section 2.1.1. We 
stress again that all matrix ensembles described in this section are isotropic, 
which implies that the ordering of the matrices is irrelevant for all statistical 
properties of the eigenvalues, see [27]. 


2.2. Singular Values 

In this section, we seek the statistical properties of the singular values rather 
than of the complex eigenvalues for some random product matrix. Explicitly, 
we consider products of Ginibre, inverse Ginibre, and truncated unitary 
matrices. It turns out that all these ensembles are polynomial ensembles [34, 
47[ which are special types of biorthogonal ensembles [45, 46]. For this 
reason we first recall a few general properties of polynomial ensembles. Let 
Xn {n = 1,...,A^) be a set of positive variables. We are interested in a 
generic jpdf of the following form 




p./3—2 

l<i<j<iV 


n ixj — Xi) det 

^ ^ ' l<Lk<N 


wl_l{xi)] , (2.57) 


Q _2 

where {w^~ (x)} is a collection of weight functions on the positive half-line. 
The k-point correlation functions are defined by [62], 


Ruxi,.... 


Xk) = 


m 


{N-ky. 


N 

n 

n=k-\-l ' 


dx. 




(2.58) 


It follows from the biorthogonal structure of the jpdf (2.57) that the fc-point 
correlation functions including the jpdf form a determinantal point process 




with kernel 


^0=2/ 


det '^{xi,Xj) , 

(2.59) 

N-l ^ 

T-Pn{x)'lpn{y)- 

n=0 

(2.60) 


The functions Pn{x) and 'ijjn{x) must satisfy the biorthogonality relation 



dx pk{x)'il;i{x) = 


hkdk,i- 


(2.61) 
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Due to the structure of the jpdf (2.57) each Pn{x) is a monic polynomial. 
For the orthogonal functions 'ipn [x) we require that their linear span given by 
span{'^o(3;), ■ ■ ■ , 'tpN-i{x)"\ agrees with that of the weight functions Wn {x), 
span{rcQ~^(x),..., The biorthogonal functions as well as the 

squared norms hn are uniquely defined. Finally, the normalisation constant 
in the jpdf (2.57) can be written as the product following [62] 

Af-l 

= iV! n K. (2.62) 

n=0 

Hence the normalisation is completely determined by the squared norms. 


2.2.1. Ginibre matrices 

Now, we are ready to discuss products of rectangular Ginibre matrices. We 
keep the notation from section 2.1.1 and consider a product matrix H^ 
given by (2.1) where each Xj is an [N + ^'j) x [N + z/j_i) matrix distributed 
according to (2.2). A derivation of the jpdf for the singular values of the 
product matrix in the /3 = 2 case was presented in [22, 23] and explicitly 
uses the Itzykson-Zuber integration formula [70]. The absence of similar 
integration formulae for /? = 1 and /3 = 4 have, so far, restricted explicit 
calculations to /3 = 2. For some asymptotic quantities this difficulty can be 
circumvented using supersymmetric techniques [53]. 

Let Xj (j = 1,..., N) denote the squared singular values of the product 
matrix Hm, then the jpdf for the singular values of the product of complex 
Ginibre matrices is given by (2.57) with weight functions [22, 23] 


w: 


P-2/ N _ 

v,k ~ '^0,M 


ui,..., + k 


x 


(2.63) 


The biorthogonal functions corresponding to these weights were hrst ob¬ 
tained in [22, 23]. It was shown that the biorthogonal functions and the 
squared norms are given by 


Pnix) 

'<Pn{x) 


- (-l)"+fen! r T{u^ + n + l) 
(-l)^+^n!r(z/M + n + l) 


(2.64) 


E 

k=0 


(n — k)\ kl T{i'm + k + 1) 

^M,0 

I/I,.. + k 


(2.65) 









20 


M 

hn = n\Y\_^{^m + n + l). ( 2 . 66 ) 

m=l 


It is immediately seen from this representation that the biorthogonal func¬ 
tions are monic, but for further calculation it is often useful to rewrite them 
in terms of special functions. The polynomial Pn{x) can be written either 
as a hypergeometric function or as a Meijer G-function, while 'ipnix) can be 
written as a single Meijer G-function. Explicitly, we have [22, 23] 


Pn{x) 


ipn{x) 


ho 


-hr,.G 


0,1 


—n 

M + 1) • • •) + 1 

n -|- 1 




X 


G 


M+1,0 


—n 




X 


(2.67) 

( 2 . 68 ) 


This determines the normalisation via (2.62) and all correlations via (2.59) 
with kernel (2.60). Furthermore, the explicit formulation of the biorthogonal 
functions (2.67) and (2.68) allows a double contour integral representation 
of the kernel [40] , 


KfGx,y) 


(27ri) 




/C 

■y0,1 


^uy-v-i y[u-n + 1 ) r{u + 1 ) 


M 


du (b dv- , . . , 1 1 . 

s u-v r(r; - iV + 1) r(r; + 1) r(r;-I-i/m + 1) 


n 


r(rt -|- I'm + 1) 


N 






-N 

0 , 1 / 1 ,..., 


uy 


(2.69) 


where C is a straight line from — ^ — ioo to — | + ioo while S encloses 
1,... ,iV in the positive direction without any intersection with C. This is 
seen by writing the biorthogonal functions (2.67) and (2.68) as their integral 
representations, see (2.10), and inserting them into the expression for ker¬ 
nel (2.60); the sum can be performed and the first line of (2.69) is obtained. 
In the second line we have displayed a second real integral representation 
from [40]. It is already very reminiscent of that large-asymptotic result 
as we will see in section 3.2.1. 

It was pointed out in [40] that the polynomials (2.64) satisfy a stricter 
condition than biorthogonality; they are multiple orthogonal polynomials of 
type II with respect to the weights w^^^{x),... :W^~M-iix) (see [71, 72] for 
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a discussion of multiple orthogonal polynomials). This means that 

) 

dz x^Pn{x)wl~^ {x) = 0 , (2.70) 


f 


for /c = 0,..., M — 1 and £ = 0,..., \— 1, where [x] denote the ceiling 
function. This multiple orthogonality may be used to establish M + 2 term 
recurrence relations for the biorthogonal functions: 


fOO 

XPnjx) = Pn+ljx) + am,nPn-m{x), flm.n = / dx X Pn{x) 
m=0 

Tpnix) 'lpn-l{x) ^ Ipn+mix 


m=0 

M 

hn-1 h. 

m=0 


n+m 


■ 5 bm,n — 


'4^n—m (^) 
hn—m 

'(Pn{x) 


I dxxpn+m{x, 

0 


(2.71) 


These recurrence coefficients were explicitly calculated in [40] but will not 
be repeated here. 


2.2.2. Ginibre and inverse Ginibre matrices 

Like for the complex eigenvalues we will also consider a mixed product of 
Ginibre and inverse Ginibre matrices; here we follow [48]. We seek the jpdf 
for the squared singular values of the matrix 

(2.72) 

where Xi and Yj are distributed according to the induced density (2.6) with 
indices i^i and pj, respectively. Recall that any rectangular structure of the 
product matrix can be incorporated by choosing Vi (i = 1,... ,M) and pj 
{j = 1,..., L) to be positive integers. The jpdf for the mixed produet of Gini¬ 
bre and inverse Ginibre matriees is given by (2.57) with weight function [48] 


w 


0=2 


{x) = G 


M,L 

L,M 


f-N - pi,..., -N - PL 

Vi,.. .,UM-i,nM + k 



(2.73) 


This choice of weight functions requires that M > 1 but allows L = 0. 
However, there exists an alternative choice which requires that L > 1 and 
allows M = 0, leading to the same correlation functions. 
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Similar to the pure Ginibre case, the biorthogonal functions may be 
expressed neatly in terms of special functions [48], 


Pn{x) 


1pn{x) 




f-n, I - N - - N - 

\ + 1, . . . , I'M + 1 



h 


/n + l,A^ + /ii,...,A^ + /ii 



*^L+1,M+1 


f-N - /ri,..., -N - fiL, -n 

0, z^i,... ,z/M 



(2.74) 

(2.75) 


M L 

hn = n\Y\ r(n + t'm + 1) r(iV + ^n-n). 

m=l £=1 


(2.76) 


If L = 0 then the biorthogonal functions reduce to the Ginibre case and 
they satisfy the M + 2 term recurrence relations (2.71). However, if L > 1 
then the multiple orthogonality relations (2.70) are no longer valid and no 
recurrence relations are known to date. 

As in the Ginibre case, the normalisation constant and the /c-point cor¬ 
relations are determined by the biorthogonal functions (2.74) and (2.75) 
together with the squared norms (2.76) via (2.62) and (2.59), respectively. 
Also a double contour integral representation is possible. 




{x,y) 


1 ^ r(u - iv +1) r(u +1) 

/s u-v r(u-iv + i)r(u + i) 

TT r(u + r'm + 1) TT r(A^ + flrn-u) 

n^r{v + Um + l)l\nN + „^-v)- ^ 


The contour S is defined as for products of Ginibre matrices and a similar 
single real integral representation as in (2.69) exists. 


2.2.3. Truncated unitary matrices 

We now turn to products of truncated unitary matrices. The derivation 
of the jpdf for such products requires an extension of the Itzykson-Zuber 
integral which was very recently proved in [52]. Moreover, the authors of [52] 
showed that the jpdf for the product of M truncated unitary matrices was 
given by (2.57) with weight function 


w 


/3=2 


(x) 


'^M,M 


/Ki,..., km-Ij f^M — N -\- k + 1 

Vi,... + k 



(2.78) 
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Here we use the same notation as in section 2.1.3. The biorthogonal functions 
are given by 


Pn{x) 


m+iFm 


-n, ki + 1,...,km + l 


1 ^0,M+1 


/-Ki,..., -KM,n + 1 



Ipnix) = G 


M+1,0 

M+1,M+1 


/-n, Ki, ...,km 


M 

hn = n\'^ 
m=l 


r(r'm + n + 1) 
T{Km + n + 1)’ 


(2.79) 

(2.80) 

(2.81) 


which determines the normalisation (2.62) and the correlations (2.59). Sim¬ 
ilar to the two previous examples, the kernel can be written as a double 
contour integral, 




1 j; ^ x^y-^-^r{u-N+l)r{u + l) 

(^y_i_,o, "'“Js n-v r(u-iV + l)r(u + l) 

TT r(tt -I- r'm -I- 1) r(u + Km + F) Qn\ 
^\T{v + Vm + l)T{u +Km + iy 

m=l ' ^ ' 


where the contour S is chosen such that it encloses 1,..., without inter¬ 
secting the other contour. For the real integral representation corresponding 
to (2.69) we refer to [52]. 


2.2.4. Mixed products 

As a final remark, we mention that more complicated products constructed 
from Ginibre, inverse Ginibre, truncated unitary and inverse truncated uni¬ 
tary matrices can be studied without introducing new techniques. Such a 
product of Ginibre matrices times a single truncated unitary matrix has 
previously been studied in [34]. When mixed products are considered it is 
important to note that all matrix ensembles described in this section are 
isotropic, which implies that the ordering of the matrices is irrelevant for 
the statistical properties of the singular values, see [27]. 


3. Local universality for large matrix dimensions 

In this section we will discuss the universal limits for products of a finite 
number M of matrices as the matrix dimension N tends to infinity. Global 
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spectra for both the complex eigenvalues and the singular values of product 
matrices have been studied intensively in the literature; in particular using 
techniques from free probability, see [13, 73, 74, 75, 12, 16, 76, 14, 17, 8]. 
Many of these results predate the exactly solvable cases described in sec¬ 
tion 2. However, almost nothing was known about the local universality for 
products of random matrices until recently. One of the main benefits of the 
matrix products described in section 2 is that their exact solvability allows 
a direct study of both global and local universality. Here we will restrict our 
attention to results about local universality. In the known cases, it turns out 
that local correlations in the bulk and at the soft edges correspond to clas¬ 
sical universality results from random matrix theory; that is Ginibre-type 
correlations for the complex eigenvalues and correlations given in terms of 
the sine and Airy kernels for the singular values. At the origin and at the 
hard edge new universality classes arise. In order to make this review both 
short and concise, our main focus will be on these new universality classes. 
We stress that this choice does not mean that the results for the bulk and the 
soft edges are less interesting, neither that these results are easily obtained. 


3.1. Complex eigenvalues 

3.1.1. Origin 

We consider the local correlations at the origin for a product of complex 
(/3 = 2) Ginibre matrices first obtained for square matrices in [19]. The 
correlations are given by (2.16) with weight (2.8) and kernel (2.19). The 
local scale at the origin is obtained by keeping the number of matrices M 
as well as the parameters Vm > (^u = 1,...,M) fixed, while taking 

the matrix dimension N to infinity without further rescaling. This is a 
trivial task since the weight function (2.8) is independent of iV, while the 
kernel (2.19) simply becomes an infinite sum which can be written either as 
a hypergeometric function or a Meijer G-function [19], 


^origin 


= J™. t-tM 




{uv*y 


N^oo 

1 + l,...,nM + l 

TT 


1 


^ nm=l r(^'m + n + l) 


UV 


n"=ir{i^~ + i) 
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= - g!’^ 


TT 




— UV 


(3.1) 
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The universal k-point correlation functions at the origin therefore read 


■ ■ ■, = n , det 

i=l 




K, 


Mm 


origin 


{Zi,Zj 


with kernel (3.1) and weight 


w 


M t 


_ ^M,0 







(3.2) 


(3.3) 


For M = 1 and 1 ^ = 0, the product ensemble reduces to the classical Ginibre 
ensemble. Evaluating the special functions (3.1) and (3.3) for M = 1 and 
inserting this into (3.2), we explicitly see that 


M=l,u=0 / \ 

PoTigin-,k 


det 


1 

— exp 
vr 



2 I 1 




(3.4) 


These are the known correlation functions for a single complex Ginibre ma¬ 
trix at the origin. Because in the Ginibre ensemble the density is flat the 
origin is not special and eq.(3.4) agrees with the bulk scaling limit, see the 
next subsection 3.1.2. 

For M > 2, the correlations (3.2) satisfy a reduction relation 


• • • , V^Zk) = Porig^ki^l, ■■■, ^fc). (3.5) 


We recall that the parameters {m = 1,..., M) incorporate the rectangu¬ 
lar structure of the matrices. 

Let us now turn to a mixed product of Ginibre and inverse Ginibre 
matrices as discussed in section 2.1.2. Here, the fe-point correlation functions 
are given by (2.16) with weight (2.42) and kernel (2.17), where the squared 
norms are given by (2.44). It can be seen that with and firn fixed for 
m = 1,..., M, we have the following scaling limit for the A:-point correlation 
function: 


lim 




Zk 


N^co origin;fc j^L/2 ’ ' ' ' ’ J^L/2 


_ M,iy 
P origin ;fc 


{zi,...,Zk). (3.6) 


Here the right hand side is given by (3.2), i.e. it is the same as for the 
product without any inverse matrices L = 0. 

Finally, we will look at products of truncated unitary matrices. We 
consider the weight function (2.50) and the kernel (2.17) with the squared 
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norms given by (2.53). Let J and L be integers such that J + L = M. We 
take M and J as well as Vm (’tt = 1, • • •, M) and Kj (j = 1,..., J) to be 
hxed, while ki = N + 0(1) as N tends to infinity ioi I = 1,... ,L. The 
fc-point correlation functions (2.16) have a hard edge scaling limit given by 


M,J,u,k/ \ _ V 


hm ^ pM,J,iy,K ( \ 

TV^oo V ’ ■ ■ ■ ’ J 


=n<-''("'),<i5L 


-^origin 


(3.7) 


where the weight and kernel are given by 


<L(^) = gJS 
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7-1 I 1, ftl + 1, Kj + 1 
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uv 


origin 


vr n™=ir(z^7n + i)/n;=ir(Ac, + i) 


\J 




0, -Kl, . . . ,-Kj 


vr ■^+bA^+i\0, 


— uv 


(3.8) 


(3.9) 


This was first shown in [32] for a product with J = 0; in this case the 
correlation functions (3.7) reduce to the case with M Ginibre matrices (3.2). 
ft remains to find the local correlations for more general products. It would 
be natural to look at mixed products constructed from both Ginibre and 
truncated unitary matrices as well as their inverses. 

* * * 

A more challenging task is to find the local correlations at the origin for 
products of real (/3 = 1) and quaternionic {(5 = 4) matrices. For /3 = 4, the 
correlations at the origin were obtained for a single matrix as well as a prod¬ 
uct of two matrices in [66] and [58], respectively. The main idea presented 
in [66] and reused in [58] was to write down a set of coupled differential equa¬ 
tions for the kernel at the origin. This technique can be extended to products 
of an arbitrary number of matrices. However, the complicated structure the 
differential equations for arbitrary M has prevented an explicit evaluation 
for M > 3, so far. The calculations simplify considerably if we integrate 
out the angular dependence of the eigenvalues; this idea was explicitly used 
in [21, 27, 30]. The local radial density at the origin for products of Ginibre 
matrices was explicitly calculated in [21] and a structure closely related that 
of complex matrices was found. For /3 = 1 the correlations for M = 2 were 
found in [59], but almost nothing is known for M > 3. 
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3.1.2. Bulk and soft edge 

It was shown in [19, 32, 33] that, under proper rescaling, the local k-point 
correlation functions of the eigenvalues in the bulk for a product of complex 
Ginibre or truncated unitary matrices are given by 



det — eyi-p i -\\zi\^ — \\zj\^ + ZiZ*) , 
i<i,j<k Itt ^ ^ ^ \ 


(3.10) 


which is the same as for the standard complex Ginibre ensemble, see e.g. [62, 
77]. In order to study the local correlations in the vicinity of a (soft) edge, 
we need to choose a point located on the edge. Given such a point, zq, it 
was shown [19, 32, 33] that proper rescaling leads to 




det 


27r 


exp ( — + ZiZj) erfc 


ZgZi + ZjZo 


, (3.11) 


where erfc(x) is the complementary error function. Note that this supple¬ 
ments the universality results at strong non-Hermiticity known for the com¬ 
plex eigenvalues of non-Hermitian matrices [78, 79, 80], see also [81] for a 
recent heuristic approach. 


3.1.3. Weak non-unitarity limit 

In the case of truncated unitary (or orthogonal) matrices one may consider 
the particular limit, in which the number of truncations remains finite while 
the matrix size(s) go to infinity. Gonsequently the resulting truncated matri¬ 
ces become almost unitary, with the macroscopic density of complex eigenval¬ 
ues condensing on the unit circle. However, locally the complex eigenvalues 
may still extend inside the unit disc. For M = 1 this limit was studied first 
in [34] and named weakly non-unitary. For M > 1 and a particular choice of 
parameters this has be generalised in [32]. Namely if in subsection 2.1.3 we 
truncate all matrices in the product starting from the same size, Kj = K, 
down to square matrices of size Nj = N (uj = 0) for all j, we take the 
following large-limit with K — N = k fixed. In order obtain non-trivial 
local correlations inside the unit disc we take k points zi in the vicinity of a 
hxed point zq on the unit circle: 

Zj = 1- ^{xj -k iyj) , Xj >0 for j = 1, 2,... , A:. (3.12) 

Here we have chosen zq = 1 without loss of generality, due to the rotation 
invariance. In contrast to the previous sections here the weight (2.50) and 
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the kernel (2.17) with norms (2.53) do not converge to a limit individually. 
Only their combination as it appears in the correlation functions (2.16) has 
a limit which we give here straight away: 


lim 

N^OO 


...A-^{xk + m)) 

j X r t^m,k / , . 


= det 


{Xj +iyj, XI+ iyi) , (3.13) 


where the weak kernel is given by [32] 


XI+iyi) = 


Q{xj)Q{xi){4:XjXi)^^^ 


- 4 ) 


tt{Mk — 1)! 

Mk / -I —t 
1 — e 


(3.14) 


t 


t={xj+xi+i(yj-yi))/M 


This expression generalises the kernel at weak non-unitarity given in [34] for 
M = 1. In order to extend this result to a more general parameter setting 
with the Kj being different, the difficulty is to obtain the limiting kernel (and 
not the weight). 


3.1.4. Further limits 

In addition to the various limits in different parts of the spectrum we dis¬ 
cussed so far further asymptotic limits can be considered. For example the 
large radius limit r —)• oo of the hole probabilities eqs. (2.23) and (2.34) 
at finite and infinite N have been derived in [20] and [30] for /3 = 2 and 
4, respectively. Furthermore the infinite-A^ process can be considered di¬ 
rectly and for example overcrowding estimates can be made [20, 30]. Similar 
questions arise in the analysis of zeros of Gaussian analytic functions, cf. 
[36]. 

A related question was considered in [38], where the distribution of the 
largest radius was studied in a double scaling limit with both N and M 
going to inhnity. In the case of products of square Ginibre matrices with 
/3 = 2 a transition between the known Gumbel distribution for M = 1 [82] 
and a lognormal distribution was found. The latter relates to the study of 
the largest stability exponent which will be introduced in section 4. 


3.2. Singular values 

3.2.1. Hard edge 

In this section we return to statistical properties of the squared singular val¬ 
ues for products of random matrices. The products described in section 2.2 
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were all determinantal point processes (2.59) with a kernel (2.60) constructed 
from a set of biorthognal functions. A study of the local correlations for such 
products at the hard edge for large matrix dimensions was first undertaken 
in [40], where it was found that a new family of universal correlation kernels 
arises. We saw in section 2.2.1 that the kernel for the product of Ginibre 
matrices can be written as a double contour integral (2.69). Using that 


r(u-iv + i) 

we find the hard edge limit 




sin TTV 




1 


1 r^M,ufX y 

A i 


(27ri)2 7_i_. 


N' N 

D 

du 


dv 


M 

n 


r(u + Vm +1) 

r(u + Vrn + 1) 


(3.15) 


^Uy V 1 sinvru r(u + 1) 
u — v sin7rur(u + l) 


(3.16) 


Evaluating the integrals allows a representation in terms of special functions 




qFm 

= ds - 

Jo 

"1 


1^1 + 1, . . . , Um + 1 


— sx 


n^=i r(^™ +1) 


G, 


M,0 




sy 




-l^M 


sxj G, 


M,0 




sy 


(3.17) 


We will refer to (3.17) as the Meijer G-kernel [52]. The Hamilton equations 
for the gab probability in this limit were studied in [50]. 

For M = 1 the product consists of a single matrix, hence it reduces to the 
Wishart-Laguerre ensemble. It is well-known that the hard edge behaviour 
for the Wishart-Laguerre ensemble is described by the Bessel kernel, see 
e.g. [77]. This is explicitly incorporated in the Meijer G-kernel (3.17) by 


= (|) ^ dsJy{2^/w)Ju{2y/^) = 4 


where the Bessel kernel is defined as 


X 


^Bessel( 4 ®, 4 y), 

(3.18) 


-^Bessel (®! V) 


y/yJu{y/x)Jl{y/y) - VxJU'/x)JuiVy) 
2{x - y) 


(3.19) 
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Note that the x and y dependent prefactor to the Bessel kernel in (3.18) 
cancels out when calculating the correlation functions due to the determi- 
nantal structure. For M >2, the Meijer G-kernel (3.17) satisfies a reduction 
relation, 

lim VMK^’''{i'MX,VMy) = (3-20) 

Um—^OO IVXCIJCI. 

where the kernel on the right hand side depends on the remaining parameters 
Um {m = 1,...,M — 1). Recall that for the product of rectangular Ginibre 
matrices, the parameters Vm incorporate the rectangular structure. 

Similar results were obtained for mixed products of Ginibre and inverse 
Ginibre matrices in [48] and for products of truncated unitary matrices 
in [52]; the product of Ginibre matrices together with a single truncated 
unitary matrix was studied in [34]. The biorthognal functions for mixed 
products of Ginibre and inverse Ginibre matrices are given in section 2.2.1. 
For M and L as well as Vm {ixi = 1,..., M) and (-^ = ■ iL) fixed, we 

have [48] 


lim 

N^oo ^ 


y 




= K, 


Meijer 


ix,y), 


(3.21) 


where the kernel on the left hand side is given by (2.60). 

For products of truncated unitary matrices the biorthogonal functions as 
well as a double contour integral representation of the kernel can be found 
in section 2.2.3. Let J and L be integers such that J + L = M. We take M 
and J as well as Um [txi = 1, • • •, M) and Kj (j = 1,..., J) to be fixed, while 
Ki = N + 0(1) as N tends to inhnity for £ = 1,..., L. In this case, we hnd 


1 


lim 

N^oo 


K 


N 


y 


J\[L+1 ’ ]\JL+1 


= (3-22) 


where the kernel on the right hand side is a generalised version of the ker¬ 
nel (3.17) given by [52] 
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(3.23) 


Note that this Meijer G-kernel reduces to (3.17) for J = 0, exactly like the 
/c-point correlation functions (3.7) reduce to (3.2) for J = 0. 
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In addition to the product ensembles described above, the Meijer G- 
kernel (3.17) has appeared in the context of Cauchy multi-matrix mod¬ 
els [42, 83, 41, 43] and Muttalib-Borodin ensembles [46, 45]. The reap¬ 
pearance the of the Meijer G-kernel (3.17) in all these models suggests a 
much stronger underlying universality principle similar to known results for 
the Bessel kernel [84, 85, 86, 87, 88]. In the global regime such a link is 
provided by the equilibrium measure [89]. 


3.2.2. Bulk and soft edge 

The evaluation of the local statistics in the bulk and at the soft edge turns 
out to be much more technically demanding than the evaluation at the hard 
edge, and has only very recently been obtained for the product of rectangular 
Ginibre matrices [51]. There it was shown that classical random statistics 
is obtained under proper rescaling, i.e. the sine kernel in the bulk and the 
Airy kernel at the soft edge. Recall that the sine and Airy kernel is defined 
by (see e.g. [62]) 


^sinei^j y) 


sin7r(x — y) 
x-y 


and KAiTy{x,y) 


Ai(x) Ah(y) - Ai^(x) Ai{y) 
x-y 


(3.24) 

respectively. We emphasise that the derivation of the sine and Airy kernel 
presented in [51] crucially depends on known results for the global density 
of the squared singular values [75, 90, 8], and in particular an asymptotic 
expression in terms of elementary functions obtained in [91, 92, 93]. How¬ 
ever, the strategy for obtaining the local statistics presented in [51] is by 
no means restricted to the Ginibre case; in fact. Ref. [51] provides similar 
results for mixed products of Ginibre and inverse Ginibre matrices (see sec¬ 
tion 2.2.2) and for Ginibre matrices mixed with a single truncated unitary 
matrix (see [34]) although the proofs are only sketched. 


4. Lyapunov and stability exponents for large products 

Until now we have considered the asymptotic behaviour of matrix products 
with a finite number of factors M as the matrix dimensions tend to infinity. 
In this section we will consider the opposite situation, where the matrix 
dimensions are kept fixed as the number of factors M tends to infinity. In 
general, we are interested in the spectral properties of a product matrix 

Hm = Xm • • • -^1, (4.1) 

where each Xm {m = 1,...,M) is an x A^ random matrix, indepen¬ 
dently chosen from some ensemble. The multiplicative ergodic theorem of 
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Oseledec [94, 95] states that if the second moments of the diagonal entries 
of XmXm are hnite, then there is a well-defined limiting matrix 

hm (4.2) 

M-^oo 

with real eigenvalues (n = 1,..., N). Here, the are known as Lya¬ 
punov exponents. Let Xn (n = 1,..., N) denote the squared singular values 
of the product matrix (4.1); Oseledec’s theorem tells us that Xn ~ 
for large M, hence negative (positive) Lyapunov exponents represent expo¬ 
nential decay (growth) and therefore stability (instability). Moreover, we 
see that all squared singular values diverge exponentially whenever the Lya¬ 
punov spectrum is non-degenerate. In such cases, it is expected that the 
Lyapunov exponents become independent Gaussian random variables with 
fluctuations of order for sufficiently large M, see [96, 6, 97]. To¬ 

gether with the symmetry of the Lyapunov exponents under permutations 
this statement is equivalent to saying that the jpdf for the squared singular 
values is given by a permanental point process. 


N 

n=l 


N\ per 

l<i j'<W 



2 ^* J 


(4.3) 


with means A* and variances (Jj to be determined. The Gaussian approxima¬ 
tion is expected to be valid when the distance between Lyapunov exponents 
is large compared to the size of their fluctuations, i.e. |Aj — Aj| S> for 

all i / j. 

How much of this expectation has been realised? In [98], the mean and 
variance of the largest Lyapunov exponent has been computed for products 
of real Ginibre matrices, see (4.4) below. The same ideas were used to 
calculate all the Lyapunov exponents in [99]. In these papers probabilistic 
methods were applied, which we will briefly recall below. Further explicit 
results for the Lyapunov exponents (again based on probability theory) were 
derived much later for factors from correlated Gaussian ensembles in [54] for 
(3 = 2 and [55] for /3 = 1,4, see Aj in (4.4) for the uncorrelated case. 

It is clear from the discussion in section 2.2 that the explicit results for 
hnite-and -M should give access to the permanental structure. Indeed 
in [9] such a derivation was recently performed for the product of Ginibre 
matrices with f3 = 2. There the permanental point process (4.3) was derived 
starting from (2.57), dehning xj = and taking M to be large. The 
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following values for the variances (and means) were obtained for (3 = 2: 

A'! = ilog| + i^(^) a.,d = (4.4) 

where 'ip{x) denotes the digamma function. For further details including the 
subleading higher order cumulants for 13 = 2 we refer to [9]. 

Previously a direct access to the jpdf at finite M and N was unknown, 
and alternative methods for the evaluation of the Lyapunov exponents were 
developed much earlier [98]. These are valid for the more general class of 
isotropic matrices (including the Ginibre case) which is why we recall them 
here. The partial sum of Lyapunov exponents can be written as 




n=l 


/3 ± 1 

N-n+l - 2M 


M 

> sup log 


delA^xlx^A 

det At A ’ 


(4.5) 


where denotes the n-th largest Lyapunov exponent and the supre- 

mum is over all x fe matrices with entries in F = M, C,]H1 for f3 = 1,2,4. 
For large M the right hand side of (4.5) can be evaluated using the law of 
large numbers. Recently, it was pointed out in [100] that the method of [98] 
also can be used to extract the variances of the Lyapunov exponents. In this 
case, formula (4.5) is used to evaluate 


k 

“ {^k) ~ ^ ^ Cov(A^_^_,_p A^_J-_|_;^). 

n=l l<i<j<k 

(4.6) 

If the Lyapunov spectrum is non-degenerate then the covariances are ex¬ 
pected to decay exponentially at large M, hence they may be neglected at 
leading order in M, see [100]. Note that if we are only interested in the 
largest Lyapunov exponent, i.e. k = 1, then the second sum on the right 
hand side of (4.6) disappears which was used much earlier in [98] to find the 
variance of the largest Lyapunov exponent. 

Also mixed products of Ginibre and inverse Ginibre matrices as well 
as products of truncated unitary matrices can be evaluated using either the 
probabilistic methods from [98] or the exactly solvable models from section 2 
as in [9, 10]. The mixed product case can be constructed from (4.4) and will 
not be repeated here, see [100]. We can consider unitary matrices truncated 
from {N + k) X (N k) Haar distributed matrices to x sub-blocks (see 
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the discussion in section 2.1.3) leading to [100] 

' fi{K + n] 




(4.7) 

(4.8) 


These quantities can also be found by direct evaluation of the jpdf from 
sections 2.2.3 and 2.1.3. While the discussion given above involves solely 
square matrices, it is a straightforward task to include the parameters Vj, 
see [10]. 

Although Lyapunov exponents are the standard choice for the characteri¬ 
sation of stability, a different possibility exists. Rather than investigating the 
singular values of the product matrix (4.1), we could look at the complex 
eigenvalues Zn {n = of the product matrix (4.1). The absolute 

values of the eigenvalues are expected to grow exponentially for large M, 
although there is no equivalent of Oseledec’s theorem in this case. It was 
therefore suggested in [9] to parametrise the eigenvalues as Zn = 
with and On real. It was previously conjectured that the stability expo¬ 
nents driving the exponential growth (decay) are identical to the Lyapunov 
exponents whenever the spectrum is non-degenerate [11, 9, 10]. For the prod¬ 
uct of Ginibre matrices this was proved in [9, 10] starting from eqs. (2.7), 
(2.31) and (2.37), leading to (4.3) with identical means and variances (4.4) 
for the stability exponents for (3 = 2, A and 1, respectively. Note however 
that the corrections to this limit which were also computed in [9, 10] differ 
from the corrections for the Lyapunov exponents known for /3 = 2 only. 

It is natural to ask about the angular dependence of the complex eigen¬ 
values for M going to inhnity as well. For /3 = 2 the spectrum is trivially 
rotational invariant; for products of (5 = A Ginibre matrices it was shown 
in [10] that the eigenvalues have a sine squared repulsion from the real axis, 
due to the pairing of complex conjugate eigenvalues. The most interest¬ 
ing case is undoubtedly (3 = 1. For the product of real Ginibre matrices 
all eigenvalues becomes real when the number of matrices tends to infinity. 
This was first observed numerically in [56] and for N = 2, and later verified 
for general N by explicit calculations in [31]. It was found that for large N 
the probability that all eigenvalues are real can be estimated by [31] 


Prob[Vj : Zj G M] ~ 



r(l/M-f i)r(i/2)\^/'^ 
r(l/M + l/ 2 ) J 


Ar2 


(4.9) 


which tends to unity as M goes to inhnity. This was conhrmed in [10]. 
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5. Open Problems 

Several open problems can be easily identified. It would be desirable to 
allow for a wider choice of ensembles of random matrices that can be multi¬ 
plied, while maintaining the exact solvability for hnite products of M factors 
at hnite matrix size ~ N. Apart from allowing for classical Hermitian en¬ 
sembles one generalisation relevant for universality is to include invariant 
non-Gaussian ensembles. As in the truncated unitary ensemble this imme¬ 
diately drops the independence amongst the matrix elements within each 
factor. On the other hand one could also drop the independence among 
different factors in the product. Furthermore, new limiting kernels may be 
accessible from the hnite N and M results, such as by taking different double 
scaling limits. 

Taking factors from the real Ginibre ensemble or studying truncated 
orthogonal matrices still presents a challenge. But even for the truncated 
unitary ensemble only a very limited parameter range of the exact solution 
at hnite N and M has been explored in the weak non-unitarity limit. This 
and the orthogonal ensemble may provide a rich study ground, being po¬ 
tentially relevant in systems with absorption or with sources of decoherence. 
Generally speaking the plethora of new mathematical results reviewed here 
has not yet been fully exploited in applications. Especially the new micro¬ 
scopic classes at the origin labelled by M will probably require to go beyond 
standard areas such as combinatorics or telecommunications. The initial 
motivation to study toy models of dynamical systems in the beginning of 
the sixties may be very worthwhile revisiting under this new light of insights. 
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